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ABSTRACT 


Procedures for obtaining position from surface gravity obser- 
vations are reviewed and their relevance assessed in the context 
of the application of modern geodetic techniques to programs of 
Earth and ocean physics. Solutions based on the use of surface 
layer techniques, the discrete value approach, and the develop- 
ment from Green's theorem are stated in summary, the latter 
being extended to order e 3 in the height anomaly. 

The representation of the surface gravity field which is required 
in order that this accuracy may be achieved is discussed. In- 
terim techniques which could be used in the absence of such a 
representation are also outlined. 

The role which can be played by the determination of changes 
in observed gravity to a few microgal, in the definition of geo- 
detic reference systems for long period studies in Earth physics, 
is discussed and the consequences of changes of zero degree 
summarized. The possible use of these techniques in future 
geodetic practice is also assessed. 
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POSITION FROM GRAVITY 


1. INTRODUCTION 
1.1 PREAMBLE 

At first glance, it would appear that geodesists of today should be grateful for 
the activities of exploration geophysicists which have made significant contribu- 
tion to all the well known gravity data banks available at present. After all, it 
is only in the past two decades that the course of events has looked favorably on 
the collection of gravity information with solely geodetic objectives in view. All 
the available surface gravity information has still not been able to provide 
meaningful definitions of position on its own at the time of writing. The tech- 
niques of satellite geodesy revolutionized practical determinations related to 
position from a consideration of the Earth's gravity field and there is little argu- 
ment that these methods give relevant position- related parameters with a pre- 
cision of about 2-4%, the higher precision estimate being obtained on the use of 
combination techniques incorporating surface gravity information with satellite 
data. 

Many factors have changed since the advent of the satellite era less than two 
decades ago. A technical advance of importance is the development of the sur- 
face ship gravimeter which provides a means of defining the surface gravity 
field in ocean areas with a resolution which if used advantageously, can be shown 
to be adequate for all present day requirements in Earth and ocean physics. A 
second development of great significance is the sensational improvement achieved 
in the precision with which absolute determinations of gravity based on inter- 
ferometric techniques, can be made (e.g., Cook 1965; Faller 1965). The perma- 
nent installation maintained by Bureau International des Poids et Mesures (BIPM) 
at Sevres, France, has been achieving a measuring precision of ±3 /xgal as a 
matter of routine for some years now (Sakuma 1971), improving the resolution 
of g from 2 parts in 10 6 to about 3 parts in 10 9 in less than 10 years. To this 
should be added the capability which has been available for some time, and en- 
ables the measurement of differences in gravity with an accuracy of at least 
1 part in 10 s on land without any appreciable measurement time. 

These significant improvements in metrology pose a series of interesting 
problems which must be dealt with before the maximum geodetic information 
can be obtained from the use of surface gravity measurements. In the first in- 
stance, it becomes necessary to review the implications of adopting a rigid body 
model for the Earth as the basis for computations of position from surface 
gravity data. A further effect to be considered are changes in the Earth space 
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location of the rotation vector and their influence on the determination of position 
from gravity. Short period mass changes smaller in magnitude than Earth tide 
effects, and possibly more difficult to model, may also have to be considered. 

Into this category fall changes in atmospheric circulation patterns from some 
model, variations in the water table and similar phenomena. Over a longer time 
scale, it is necessary to consider the implications of a possible secular change 
in the gravitational constant G. 

As most of these effects are 7-8 orders of magnitude smaller than that of g, it 
has been accepted practice to convert observed gravity g to the gravity anomaly 
Ag by differencing g from the value of normal gravity y for a model of the Earth, 
afforded by the value of GM, the rate of rotation a> of the Earth, assumed to be 
constant, together with the equatorial radius a and flattening f of an ellipsoid of 
revolution which "best fits" the geoid. No allowance is made for the possibility 
of variations with time, in any of the parameters defining the system of refer- 
ence. This is not inconsistent with the concept of determinations relevant to a 
certain epoch, provided 

(i) the observations used are all made during the epoch considered; and 
(ii) the accuracy sought is less than 1 part in 10 7 in g. 

The order of magnitude of gravitational deviations from a solid Earth model are 
smaller than o { 10 - 6 g} . The largest effect is the diurnal Earth tide variation 
with magnitude o (10- 7 g }. It has not been considered necessary at the present 
time, to recommend the adoption of a systematic procedure for modeling and 
removing the effect of Earth tides from observed gravity except when establishing 
gravity standardization networks, in view of the limited accuracy of elevation data 
used in computing the gravity anomaly. This would call for the acceptance of a 
universally acceptable model for Earth tides, which would be used as a matter 
of routine to correct observed gravity prior to use in geodetic computations. 

Such corrections are only necessary at fundamental gravity stations at which 
determinations are made with the highest possible accuracy for either the defi- 
nition of the global gravity standardization network (Mather 1973, p. 68) or when 
attempting to locate changes in the position of the Earth's center of mass (Geo- 
center) with time (Mather 1972, p. 13). The need for applying such corrections 
at other stations will depend on the extent of gravity coverage available globally 
and whether the elevation datums have been unified at the 50 cm level. 

Current practice accepts the validity of each individual nation's elevation datum 
as well as its gravity datum. The continuance of such a practice is unwise if 
systematic errors at the 50 cm level are not to occur in the final results. The 
most taxing goal in the definition of position from surface gravity, is the deter- 
mination of the geoid in ocean areas to the highest possible accuracy, in order 
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that such results could be used with satellite altimeter data to recover the sea 
surface topography. On present trends, it would appear that 5-10 cm accuracy 
is desirable in the geoid determinations for a meaningful evaluation of the sea 
surface topography. The determination of the latter is vital for the study of 
tides in open oceans, ocean circulation and stationary characteristics of sea 
surface topography. 

It is in this context that the use of gravimetric techniques in the determination 
of geodetic position should be reviewed. The present development covers 

(a) the basic principles underlying the determination of position from 
gravity; 

(b) a review of some of the methods suggested to the present time, for 
solving the boundary value problem in physical geodesy; 

(c) techniques for the preparation of data sets for this task; and 

(d) the geodetic interpretation of such solutions in Earth space. 

In all sections, an assessment is made of the requirements which will have to 
be met in order that the independent evaluation of selected geodetic characteris 
tics available from surface gravity determinations, can be used in the resolution 
of some possible ambiguities from other methods when applied to high precision 
studies in Earth and ocean physics. 


1.2 A GUIDE TO NOTATION 
1.2.1 Recurring Symbols 

a = equatorial radius of the ellipsoid of reference 

d = separation vector between equivalent points P on the Earth's surface 
and Q on the telluroid 

dz = increment in orthometric elevation 

do- = element of surface area on unit sphere 

e = eccentricity of the meridian ellipse; e 2 = 2f - f 2 

F(i p) = f (0) sin 0 

f = flattening of meridian ellipse 

f{0) = Stokes function = cosec 1/2 1 // +1-5 cos 0 - 6 sin 1/20 - 
3 cos 0 [log (sin 1/2 0 (1 + sin l/20)jj 
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g = 
h = 



M = 


M {X> = 

m = 


m' = 
N = 
N = 



R = 
R = 


R b =* 

R m = 

r = 


r = 

= 

S = 
U = 

U 0 = 

v d - 

w = 
w„ = 


X. 

1 

X. 


observed gravity at the surface of the Earth 
ellipsoidal elevation 
height anomaly 

normal height, defined prior to Equation 6 
mass of the Earth, including the atmosphere 
global mean value of X 

aw 2 / y e 

a 3 <^> 2 / GM = m + o{f 2 } 

elevation of geoid above ellipsoid 

unit vector normal to the surface of the Earth 

Free Air Geoid; the Stokesian contribution to h d 

Indirect effect to free air geoid; non-Stokesian contribution to h d 

geocentric distance 

radius of sphere containing all topography (Brillouin sphere) 

radius of sphere which is internal to the Earth's surface (Bjerhammar 
sphere) 

mean radius of the Earth 

distance between the point of computation P and the element of surface 
area dS 

2 R sin 1/2 0 

2 R sin 1/2 0 

surface of the Earth 

spheropotential due to the system of reference 
spheropotential on the surface of the reference ellipsoid 
disturbing potential 
geopotential 
potential of the geoid 

geocentric rectangular Cartesian co-ordinate system Xj X 2 X 3 

local rectangular Cartesian co-ordinate system x t x 2 x 3 with the x 3 
axis along the local normal, the x x x 2 plane defining the local horizon, 
with axes oriented north, east respectively. 


4 



a = azimuth 

fi = ground slope; subscripts 1 and 2 refer to components north and east 

7 = normal gravity due to the reference system; subscript 0 refer to values 
on the reference ellipsoid; subscript e refers to equatorial value. 

A g = gravity anomaly at the surface of the Earth, defined by Equation 10. 

AW = geopotential difference with respect to the geoid 

S g = gravity disturbance 

A. = longitude, positive east 

f = components of deflection of the vertical; subscripts 1 and 2 refer to 
those in north and east directions respectively 

£ = deflection of the vertical, positive if outward vertical lies north, east 
of normal. 

$ = density of surface layer, except in section 3.4. 

4> = latitude, positive north; subscripts c, a and g refer to geocentric, 
astronomically determined and geodetic latitudes respectively 

0 = angular distance at geocenter between the point of computation P and 
element of surface area dS. 

a; = angular velocity of rotation of the Earth 





1.2.2 Conventions 

a = b + o {b 2 } = terms whose order of magnitude is equal to or less 

than b 2 are neglected (b < 1 ) 

x a y a = x i y i + x 2 y 2 

x i = a ijkj “ x i “ a ii fc> x + a i2 b 2 + there being as many 

equations as possible values of i 

a = c = a is approximately equal to c 
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2. BASIC PRINCIPLES 


2.1 THE SYSTEM OF REFERENCE 

The determination of Earth space position from gravity observations at the 
surface of the Earth, as pointed out in section 1, is implied from deviations of 
observed gravity from those at an "equivalent" point on some Earth model, 
whose parameters are completely defined. Current geodetic practice (LAG 1970, 
p. 12) specifies a rigid body model by the following parameters. 

(a) The value ^ (= G M) where G is the gravitational constant and M the mass 
of the Earth. 

(b) The constant rate of rotation a> of the rigid Earth model. 

(c) The equatorial radius a of an ellipsoid of revolution which presumably 
is one of best fit to the geoid. 

(d) The dynamic form factor J 2 which is equivalent to a value of a flattening 
f for the reference ellipsoid. 

It is conventional to choose a reference ellipsoid whose equatorial radius is such 
that it has the same volume as the geoid. This is not a necessary condition if 
zero degree effects are taken into account when formulating a solution for the 
boundary value problem. What is more important in solutions which aspire to 
accuracies greater than the order of the flattening (i.e. , 30 cm in the height 
anomaly), is that the ellipsoid lies everywhere within the physical surface of the 
Earth. This enables the use of Laplace's equation in the representation of the 
appropriate disturbing potential without approximation. 

The adoption of such a procedure without an equivalent adjustment in fj. could 
cause larger values of the gravity anomaly which would in turn, call for greater 
caution in developing computer algorithms for numerical work. 

It can be stated without being contentious that the value adopted for a has to be 
based on some determination of the scale of Earth space. This would be provided 
by either the measurement of long arcs at the surface of the Earth by classical 
techniques (e.g. , the Pageos baselines) or else by laser ranges to either satellites 
or the moon. All determinations of scale are therefore based on the velocity of 
light. 

The value of the flattening f of the reference ellipsoid is best deduced from the 
second degree zonal harmonic obtained from the secular variations in the right 
ascension Q of the node and the argument co of perigee of near Earth satellites. 
The precision claimed at the present time (e.g., Lerch et al 1972, p. 27) for this 
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harmonic is 1 part in 10 s , the required relation being (e.g. , Mather 1971a, 
P. 85) 


C 


20 




+ o{f 3 } , 


(1) 


where 


a 3 ^ 2 

GM 


co being the angular velocity of rotation of the Earth. 

The exact relation between the observed secular variations fi 20 and C 2o is (e.g., 
ibid, p. 151) 


O - 3( GM ) 1/2 _*!_ . r 

- - cos i C 


20 


2(l-e s 2 ) 2 aj/ 2 


20 > 


(2) 


a s being the equatorial radius of the satellite orbit, e s its eccentricity and i its 
inclination. 


The change df in f due to changes da, d.w and d(GM) in a, co and GM are given 
by 


df 




d co 



1_ d (GM) 
2 GM 


(3 m' + C 20 ) + o{f 2 df}. 


( 3 ) 


The ratio df/ f is therefore of the same order of magnitude as d(GM)/ GM for 
a specified value of a as the ratio dco/co is at least an order smaller if these 
ratios reflect the precision with which each of the quantities are determined. 

It is all important that the rotational characteristics assigned to the reference 
model are exactly equivalent to those influencing gravity as measured at the 
Earth's surface. This is implicit in deriving Equation 60 from Equations 58 and 
59. The rotation vector in Earth space is not fixed and hence deviations from 
the rigid body model adopted for any system of reference (e.g., IAG 1970, p. 24). 
The rate of rotation co is subject to secular variations due to the effects of tidal 
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friction, and certain short period effects which, in theory, have to be accounted 
for when reducing observed gravity to a rigid body equivalent of the Earth, are 
not of practical consequence as they have a magnitude of less than 1 part in 10 9 
in g. A second factor is the change in position of the instantaneous axis of rota- 
tion with respect to the Earth’s crust. The total contribution to observed gravity 
of the rotation is the appropriate resolute of 


g r - P" 2 ( 4 ) 

directed away from the axis of rotation and perpendicular to it. The changes in 
g r due to c hang es dp in p, which is the distance of the point at which gravity is 
measured from the axis of rotation, and dw in w is given by 


g r 



dp 

P 


+ o{ lO- 



tS) 


The effect of the ratio dw/o; will be less than 1 p gal on observed gravity for a 
10 msec variation in the length of day and hence this term is not of significance 
in the reduction of observed gravity, if the latter were restricted to some epoch 
of observation. The effect of polar motion is o { 1 /xgal) (Bursa 1971). 

More significant short period changes in observed gravity have been reported by 
Sakuma (1971) after modeling the effect of Earth tides. It should be noted that a 
1% change in the local atmospheric density is of o {10 /xgal } while quasi- stationary 
changes in the local geological formations, e.g. , in the local water table, could 
cause gravitational effects of this same magnitude. It is therefore important that 
both the atmosphere and the local geology be modeled in the vicinity of those 
gravity stations at which g is to be re-measured at intervals of time with the 
highest possible precision. 

The dominant gravitational variation with time is that due to Earth tides, with 
magnitude of o {10 2 pga \ } and Earth models for this effect are well known in the 
literature (e.g., Melchoir 1966). It is important that an unambiguous Earth tide 
model at the 10 p gal level be uniformly adopted when specifying gravity values 
at stations comprising the global gravity standardization network described in 
section 4. 

The final parameter defining the reference system is p(- GM). The commonly 
accepted values of GM are all based on the analysis of interplanetary space 
probes. The technique used can be briefly summarized as follows (Esposito 1972). 



Doppler data from interplanetary spaceprobes is analyzed using numerical 
integration procedures for the determination of the motion of the probe with 
reference to a geocentric inertial co-ordinate system. Perturbations due to the 
Earth's departure from a sphere, solar radiation pressure, planetary and lunar 
gravitational effects and spacecraft attitude control forces are modeled when 
effecting this solution which also provides revised estimates of tracking station 
co-ordinates as a by product of the solution. The main conclusions of relevance 
to the present review are the following. Firstly, the value of G M is based on the 
velocity of light. Secondly, the potential U 0 on the surface of the reference 
ellipsoid, assumed to be an equipotential is related to the adopted values of GM, 
a, f and co by the relation (e.g., Mather 1971a, p. 83) 


U n 


GM 

a 


sin 1 e 


3 ^ 


where 


e 2 = 2f - f 2 . 


As GM has an uncertainty of 1 part in 10 6 at the present time, it follows that U 0 
will differ from the true potential of the geoid consistent with Newtonian gravita- 
tion, as scaled by the velocity of light by at least 1 part in 10 6 . If U 0 were as- 
sumed to be equal to the potential W 0 of the geoid, it would be tantamount to im- 
posing a second scale constraint when using gravitational techniques in geodesy. 
The only way out of this impasse is to attempt to estimate W 0 by using purely 
geometrical techniques and thereby make position determinations using gravity 
data consistent with the scale provided by the velocity of light. In the interim, 
it must be borne in mind that all position determinations based on gravity may 
have constant scale error of up to 1 p.p.m. on this account. 

In summary, it may be stated that 

(a) Only Earth tide effects need to be allowed for in all work except those 
determinations required for the monitoring of co-ordinate systems; 

(b) Position determinations based on gravity are liable to have a constant 
scale error of up to 1 part in 10 6 due to the uncertainty in G M. 

The following conclusions may therefore be drawn about the adoption of a rigid, 
body model of the Earth as an intermediary in the definition of position from 
gravity, with the highest possible precision as the ultimate goal. 
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(a) The only departures from rigidity which need be considered for solutions 
of the boundary value problem are the effect of Earth tides on gravity 
observations constituting gravity standardization networks. 

(b) An error in GM will give rise to ambiguities in scale with corresponding 
magnitude unless independent techniques are developed for using geo- 
metrical means for determining the potential W 0 of the geoid. For 
further discussion, see section 5.1. 


2.2 DATA REQUIREMENTS 

Observed gravity will be the results of two kinds of determinations. The first 
type will be absolute determinations with the highest precision possible while the 
second will be point values established by differential techniques based on the 
absolute determinations, with a precision which is at least one order of magnitude 
inferior. The practice adopted in gravimetric determinations is the use of gravity 
observations and a knowledge of the surface topography of the Earth to determine 
the separation vector d between "equivalent" points on the reference model, 
whose Earth space position is known, and the physical surface of the Earth, as 
illustrated in Figure 1. The separation vector can be completely defined by the 
height anomaly h d and the angles g a which are more completely defined in the 
next sub- section. 

The most exacting requirements are called for in the definition of position from 
gravity when determining the geoid for ocean physics applications, where present 
estimates of requirements call for resolution at the ±10 cm level. The equivalent 
order of magnitude is e 3 (i.e, , 5 parts in 10 4 ) which can be assessed as ±50 /x gal 
in the gravity anomaly Ag. On the basis of the discussion in the previous sub- 
section, it would be adequate to maintain a rotating rigid body model as the system 
of reference and apply the appropriate reductions to observed gravity to make the 
measurements compatible with the model. The position so defined will be unaf- 
fected by short period time variations in the Earth's gravitational field. 

The nature of the reductions necessary will depend on the purpose for which the 
gravity data is required, the accuracy with which it has been established and the 
nature of the elevation data available for its reduction. Corrections for Earth 
tides are necessary both for stations monitoring changes in the global reference 
system (Mather 1972, Sec. 3.3) and to an accuracy which is one order of magni- 
tude less for stations in the global gravity standardization network. The model 
adopted for Earth tides should therefore be capable of resolution to 1 jugal. Some 
difficulty may be experienced in removing ocean loading effects in coastal areas 
which influence the tidal correction in the second significant figure (Hendershott 
1972). 
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VERTICAL 



Figure 1. The Separation Vector d and the Reference System 


Corrections for short period variations in the stratification of the atmosphere 
and the stratigraphy in the vicinity of gravity stations monitoring the reference 
system are also necessary. This presupposes the existence of "accepted" models 
for both the atmosphere and the local stratigraphy, and the effects meteorological 
changes may have on them. While Earth tide effects should be allowed for when 
making any gravity determination, the summary in section 4 indicates that the 
effect of omitting the correction on determinations of position from gravity is 
likely to be negligible as the tidal effect has the characteristics of a random 
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measurement error. All further discussion will assume observed gravity g to 
be measured on a solid Earth, rotating with a uniform angular velocity, the 
gravitational effects of polar motion being allowed for when gravity data is used 
for the definition of geodetic reference systems. 

The formulation of relations at the surface of the Earth is based on the following 
principles. 

(a) An estimate is available of the geocentric co-ordinates of the point P at 
the surface of the Earth. In classical terms, these surface co-ordinates 
(^a> A a ) are re lated to the vertical at P by astronomical determinations, 
and can be determined at best to a factor or two better than 1 part in 10 6 
(i.e., ±6 m in position in each co-ordinate). This estimate differs from 
the true geocentric co-ordinate by amounts up to 10“ 4 radians, depending 
on the magnitude of the local deflection of the vertical. 

(b) The displacement of P above the equivalent point P Q on the ellipsoid is 
defined by the normal elevation h n , which is related to the difference in 
geopotential AW between the equipotential datum for elevations (the 
geoid) and P as obtained from levelling by the relation 

AW = - f P g dz, 

-'geoid 


g being observed gravity for the section of the line of levelling where 
the orthometric height difference is dz. The equation defining h n in 
terms of AW is the relation (e.g., Mather 1971a, p. 100) 


AW / Aw \ 2 

1 + TV ( 1 + m+ f- 2fsin 2 0) + ( +o{f 3 } 

7 o \ a 7 0 / 


where y 0 is the value of normal gravity on the reference ellipsoid and 


(6) 


y ? being the value of normal gravity at the equator, and AW is treated 
without regard to sign for points exterior to the geoid. 


AW 
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It has been shown (Mather 1973, sec. 4.3) that the data requirements for the 
determination 'of the height anomaly h d with a precision equivalent to that pos- 
sible in establishing h n are well within the capabilities of measuring and data 
sampling techniques available at the present time. 

Thus position determination from gravity in any absolute sense 

(a) requires a knowledge of astronomical co-ordinates; and 

(b) a global representation of the gravity anomaly field. 

The resolution of the information from positional astronomy at the present time 
will have to improve by a factor of 50 before the horizontal determinations are 
of adequate accuracy for the complete determination of position by this method 
alone. Such a determination will also require the determination of deflections of 
the vertical to o (10 -4 <f a } . 

It can therefore be concluded that the determination of geocentric position from 
positional astronomy and surface gravity to accuracies much in excess of 1 part 
in 10 may not be a practical possibility in the foreseeable future. The determi- 
nation of the height anomaly, on the other hand, will remain a problem of funda- 
mental interest as it forms an integral part in the definition of sea surface 
topography from space. The ensuing development will continue to deal with the 
complete development necessary for the determination of position from gravity, 
but only in outline. More detailed review will be confined to the techniques for 
determining the height anomaly. 


2.3 BASIC RELATIONS 

The formulation of the Molodenskii problem (Heiskanen and Moritz 1967, p. 291) 
can be treated as one which seeks the determination of the separation vector d 
between "equivalent” points P (<*> g , W p = W 0 + AW) on the Earth's surface 
and Q(0 a » k a , Uq - U 0 + AW) on the associated spherop U = U 0 of the refer- 
ence system, as illustrated in Figure 1. If the subscripts a refer to values de- 
termined astronomically determined at P, the separation vector is given by 


d “ R a£a 5 + h d 3> (7) 

where R a are the meridian and prime vertical radii of curvature of the associated 
spherop, a are unit vectors oriented along the tangent plane to the spherop at Q 
in the meridian and prime vertical respectively, while 3 is the unit vector along 
the outward normal at Q. The subscripts g refer to the surface co-ordinates of 
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I 


the point P' in Figure 1 on the associated spherop U = U Q whose normal passes 
through P. 

The locus of the point Q is called the Telluroid and mirrors the physical surface 
of the solid Earth and the oceans to order f 2 . The system of reference is based 
on the family of spherops (U = U 0 + AW) exterior to the reference ellipsoid, 
defined by the geometrical parameters a and f , and the gravitational characteris- 
tic /j. (= GM), with the constraint that the surface of the reference ellipsoid is 
the equipotential surface U = U 0 . As pointed out in section 2.1, there is no 
necessity for the ellipsoid to be forced to have the same volume as the geoid 
(W = W Q ) if terms of zero degree were retained in the solution. In such circum- 
stances, it is easy to show that (e.g. , Mather 1973, p. 14) the height anomaly 
(h d in Figure 1) is given by 


h d = 7 [ v d " -M + o{10- 3 m}, (8) 


where V d is the disturbing potential at P given by 


V 


dp 



( 9 ) 


The quantity W p is defined by the value W 0 of the geopotential on the equipotential 
surface used as the datum for geodetic levelling, and the observed difference of 
geopotential AW between this surface and P, by the relation 

W = W n + AW. 

P U 


The datum in use at the present time is that afforded by mean sea level derived 
from tide gauge readings over periods in excess of one year. As solutions of the 
geodetic boundary value problem require definitions which are applicable globally, 
it is essential that all regional definitions of mean sea level are correlated on a 
world wide basis to a common epoch in the first instance before the differences 
in geopotential (AW) can be considered to be referred to an equipotential surface 
of the Earth's gravitational field. It has been estimated that systematic errors 
of o {±10 cm) could result in solutions of the boundary value problem if errors 
on this account were of o {±30 cm} and each datum covered o {10 6 km 2 } (Mather 
1973, p. 68). 

A second problem of consequence and of which little is known at the present time, 
are the quasi-stationary departures of the sea surface from the equipotential 
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surface defined by the results of geodetic levelling. The phenomenon, known as 
stationary sea surface topography, has been reported along coastlines in many 
parts of the world. A summary of results is given by Hamon and Greig (1973), 
indicating magnitudes of 30-50 cm being commonplace, with one reported as 
large as 1.7 m. This effect is discussed further in section 4. 

The gravity anomaly Ag at the surface of the Earth is defined by 


A e - g p " V , 


( 10 ) 


where g p is the value of observed gravity at P, corrected for departures of the 
Earth from the rigid body model, as described in section 2.1, and y pl is normal 
gravity due to the reference system at P' in Figure 1. y p , is obtained from the 
equivalent value y 0 of normal gravity on the reference ellipsoid, given by the 
commonly used relations of the type (e.g., Heiskanen and Moritz, p. 79) 


y 0 - y e jl + /S sin 2 <£ g + P 2 sin 2 2 4> g + o{f 3 } 


( 11 ) 


where y e is equatorial gravity defined by the values adopted for a, f, GM and co 
(e.g., Mather 1971, p. 87; LAG 1970, p. 48), /3 = o {f} and fi 2 = o {f 2 }, using the 
relationship (e.g., Mather 1971a, p. 101) 


V 


7o - 2 


AW 


1 + f + m - 2f sin 2 4> 


1 AW fr2X 

o ~ + O { f 2 ) 

2 ay -j 


(12) 


AW having the same significance as in Equation 6. 

Possible sources of systematic error in the computed value of y and hence Ag 
arise in the definition of <t> and AW. While the effect of errors in the latter have 

' g 

already been described, should be defined to ±0.4 arcsec if y pl is not to have 
an error of ±10 /x gal. It is therefore important to use any of the global solutions 
available at the present time for the definition of geocentric position, to evaluate 
geocentric orientation parameters for each of the regional geodetic datums 
(Mather 1973, p. 16) before computing gravity anomalies for high precision de- 
terminations rather than use values referred to regional geodetic datums. 

The equation described as the fundamental equation in physical geodesy (Heiskanen 
and Moritz 1967, p. 86), defines the relationship between the disturbing potential 
V d and the gravity anomaly Ag as (Mather 1973, p. 18) 
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(13) 
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3 y 

3 h 


d ( + 2 g ^ + ° 


where the terms within the bracket take into account effects smaller than o (f Ag} , 
i being the deflection of the vertical at the point considered. 

The philosophy underlying Equation 13 is the contention that the geocentric posi- 
tion of P is not known, though estimates adequate for the linearization of the 
quantities involved, are available. Circumstances may well arise in the future 
where accurate horizontal and vertical surveys may be available and the principal 
practical role of techniques in physical geodesy is the determination of the geoid 
in ocean areas for studies of sea surface topography. In such a situation, it is 
envisaged that all the land masses are linked to a geocentric system of reference 
using laser ranging methods and/or VLBI, giving at least one fundamental station 
on each geodetic datum. Horizontal survey methods together with geodetic and 
astrogeodetic levelling will provide data for completely defining geocentric posi- 
tion of points on any regional network which includes at least one fundamental 
station, with an accuracy of 1 part in 10 6 . As surface ship locations can be rou- 
tinely determined to within one order of magnitude greater, it is of relevance to 
examine the gravity disturbance Sg (e.g., Hotine 1969, p. 312) given by 


Sg 


g P " 


B V d 

3 h + 


1 

2 


g i 2 + o {1 /zgal} , 


(14) 


where the uncertainties in defining the position of P can be estimated as ±0.2 arc- 
sec in horizontal position and ±2 m in normal displacement, if the astrogeodetic 
levelling is based on an adequate distribution of stations. The effect of errors 
dug to the first source on Sg are o { 1 ,ugal } while that of those due to the second 
are o (5 x 10 2 /zgal }. Thus the gravity disturbance, whose order of magnitude is 
not significantly different from that of the gravity anomaly, is likely to have errors 
of o (5 x 10 2 ^gal} which are probably correlated with wavelengths in excess of 
1000 km (e.g., see Mather, Barlow and Fryer 1971, Fig. 4.2) unless radically new 
techniques are available for determining 

either geocentric position at each gravity station such that the radial com- 
ponent is resolved with systematic biases of wavelengths longer 
than 1000 km held to below the 20-30 cm level; 

£r the contribution of astro-geodetic levelling with the same resolution 
as geodetic levelling. 
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The projection of present day techniques does not lead to the conclusion that 
there would be significant advantages in using the gravity disturbance Sg in 
preference to the gravity anomaly A g in formulating solutions of the boundary 
value problem. 

The separation vector d, illustrated in Figure 1, can be represented by com- 
ponents along the axes of a local Cartesian co-ordinate system x i at Q, with the 
x 3 axis oriented along the spherop normal at Q and the x t x 2 plane lying in the 
horizon at Q, as illustrated in Figure 2, in accordance with Equation 7. d is of 
importance in defining the geocentric orientation vector 0 for regional geodetic 
datums using surface gravity data (Mather 1971b, p. 62). A description of how 
such information could be used to assemble a world geodetic system linking the 
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major land masses by comparing the separation vectors as obtained from 
gravimetry and astro-geodesy is given by Mather (1971c). 

The mainstream of practical endeavors at the present time is in the determina- 
tion of the height anomaly h d . Over 90% of the power in such determinations 
comes from the "free air geoid" N f obtained by the use of free air anomalies 
(i.e., gravity anomalies to the order of the flattening) in Stokes' integral which 
is set out in Equation 16. The latter is a solution of the boundary value problem 
for a spherical Earth which is exterior to all matter and whose bounding surface 
is an equipotential (Stokes 1849). 

The deflections of the vertical are usually obtained using the principles gen- 
erally attributed to Vening Meinesz (1928). Working on a spherical reference 
system, he showed that if the separation N f between the physical and reference 
surfaces were given by Stokes' integral 

N r P = JJ fW^gdo-, (15) 


where Ag is the value of the gravity anomaly at the element of surface area da 
on unit sphere which is at an angular distance 0 from the point of computation P, 
f (0) being Stokes' function (e.g. , Heiskanen and Moritz, 1967, p. 94), then 


^ a = 


9N 


(16) 


as illustrated in Figure 3, where x a is a two-dimensional Cartesian system in 
the horizon plane at the point of computation P, with the x t axis oriented north 
and the x 2 axis east. It is not difficult to show in the case of Stokes' problem 
that 




4 77 7 J 


3 [f(0>] 

3^— cos A a Ag do-, 
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Figure 3. The Vening Meinesz Problem 


where 


A x - a and A 2 - 7^ 77 ” a > (19) 

a being the azimuth of do- from P. This follows as only 4> in the kernel of the 
integral at (15) changes as the point of computation changes from P to some 
adjacent point Q in the case of Stokes' problem. 

The required expression for the Molodenskii problem is not the same as the 
elevation h p of P appears in the kernel of the integral. As the deflection of the 
vertical at the surface of the Earth is obtained from the height anomaly h d (ibid, 
p. 312), which is given by 

h d ~ h d(^>^’ h p) ~ h d W'- a cr’ h p) >' ( 2 °) 


where a a is the azimuth of P from the element of surface area dcr, it can be 
shown that (Mather 1971c, p. 88) 


1 

"3 h d 

3 dj 

Bh d 3 a, 

K 

d 

3u. 

+ 3a cr 9 u a . 


( 21 ) 
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u a being a set of curvilinear surface co-ordinates on the reference surface, and 
h a the associated linearization parameters. For the latitude-longitude system 

u i = <t> , u 2 = (22) 


while 


N ' F and h 2 - R cos 0 (23) 

for a spherical approximation of the Earth. In the case of solutions to order e 3 , 

W " (P + h) ; h 2 = (y + h) cos 0, (24) 

P and v being equivalent to R a defined in Equation 7. 

The detailed development of solutions of the boundary value problem in this case 
use the geocentric latitude 0 C instead of the geodetic latitude 0, all parameters 
referring to quantities relating angular displacements between the geocentric 
radii to the pole, P and dcr. 

As explained above, the principal task in determining position from gravity, is 
the definition of the height anomaly h d , which is equal to the geoid height N in 
ocean areas, where Aw = 0. The next section deals in summary with some of 
the methods which have been proposed for obtaining the height anomaly. 


3- TECHNIQUES FOR THE SOLUTION 

QF THE BOUNDARY VALUE PROBLEM 

3.1 INTRODUCTION 

Attention has been confined to three techniques whose use to obtain solutions to 
the boundary value problem has been extensively reported in the literature. The 
methods considered deal with formulations of solutions to what is known as 
Molodensku's problem, at the physical surface of the Earth. It is not intended 
to formulate solutions for surfaces other than that of the Earth, e.g., the geoid, 
obtained by defining N instead of h d . Neither is any attempt being made to dis- 
cuss the merits of regularization (e.g., Molodenskii et al, p. 45 et seq.), where 
the conditions applicable to Stokes' problem are artificially created by the trans- 
fer of mass to within the geoid. The main advantage claimed for such techniques 
is a utility which is desirable when the surface gravity coverage is poor, as the 
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adoption of certain types of mass transfers enables more reliable predictions of 
gravity anomalies for the chosen model . The validity of such claims is open to 
question if the end-product of the calculations is to be a meaningful determina- 
tion of positional parameters, which can only be as good as the available data. 

The three techniques which will be covered, ostensibly do not require a knowledge 
of the stratification of matter within the Earth, defining solutions in terms of an 
"adequate" sampling of gravity at the surface of the Earth, in conjunction with a 
complete definition of the associated topography. They can be classified as 

(1) Surface layer solutions; 

(2) Solutions from data sampled at discrete points on the Earth's surface; 
and 

(3) Solutions from Green's third identity. 

It is of interest to summarize the basis of each of these methods. 


3.2 THE SURFACE LAYER TECHNIQUE 

This method, initially developed by Molodenskii (Molodenskii et al 1962, p. 118 
et seq.) was first published in 1949. Considerable material is available on the 
problems associated with the practical use of this technique by Moritz (1966; 
1970; 1972) and members of the Soviet school (e.g. , Brovar 1964; Marych 1969; 
Yeremeyev 1969; Pellinen 1972). The derivation calls for the representation of . 
the disturbing potential V d at the surface of the Earth by a surface layer of den- 
sity $ such that the former can be represented at any point P either on the 
surface of the Earth or exterior to it by the relation 

= If s 7 ds ' < 25 > 


where there is no restriction on the shape of the surface S. It can be shown that 



2 77 <1> COS 6 
p r p 



(26) 


where the subscript p refers to evaluation at P, fi being the ground slope at P. 
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The first term on the right appears due to the indeterminance at P itself. The 
inner zone in this region is treated as a disk (e.g., Heiskanen and Moritz 1967, 
p. 129), the negative sign being introduced as the outward derivative is required, 
while the attraction of the disk is toward the geocenter. The cos /3 term allows 
for slope of the surface of the disk with respect to the vertical. No approximations 
are involved in the derivation of Equation 26. 

The ensuing development, which is well documented by Heiskanen and Moritz 
(ibid, p. 300) can be summarized as follows, retaining those terms whose contri- 
butions are greater than o {fh d > . On using Equations 8, 13, 25 and 26, 


Ag = 2 77 <E> cos B ~ 

& p r p 


W 0 /3-y 
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and 


r = (R 2 + R 2 - 2R R cos 0) V2 , 

V P P r 7 


(29) 


it follows from Figure 4 that 
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0 dS + o {f Ag} 


(31) 


The solution suggested by Molodenskii (1962, p. 120) for Equation 31 is based on 
the method of successive approximations where the surface S of the Earth is 
transformed into the surface S using a parameter k which specifies the relation- 
ship between the geocentric radii R and R to equivalent points on S and S by the 
relation 


R = R m + k (R - R m ) , 


(32) 


where R m is the mean radius of the Earth, and 0 < k < 1. Thus S and S coincide 
when k = 1, while the classical Stokesian case in which no topography exists ex- 
terior to the geoid, is obtained when k = 0. This is equivalent to scaling all ele- 
vations and grades by k from h and tan ft to h and k tan j3 where 

h = k h 

and the related angle /3 is given by 

y3 = cos -1 ^[1 + k 2 tan 2 /S] -1/2 j . (33) 

Other relevant conversions are 

r = [r 0 2 + k 2 (h - h p ) 2 J 1/2 + o{f 7} , (34) 

where r 0 is the expression for the spherical case, given by 


1 

r o - 2R m sin - 2 ^- 


(35) 


Molodenskii simplifies the solution by introducing the parameter X defined by 
the equation 


r2 

X = <D sec /3, (36) 

R 2 
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which, when taken in conjunction with the relation 


dS - R 2 do- sec fi , 


( 37 ) 


where do is the element of surface area on unit sphere, enables Equation 31 to 
be written as 


Ag 



2ttX 

p 



\ - Up /dy\ 

y [dhj 

' P C 


(38) 




X do + o {f Ag} . 


It can be shown (ibid, p. 121) that if X were expanded in a power series of the 
form 

. oo 

X = £ k 1 , (39) 

i=0 

and on introducing a set of functions G the use of Equations 33, 34 and 39 in 
Equation 38 gives a system of integral equations 


W o - U 0 3 f f 

G i = 2-77 X. + 2 — J R m I — do , i = 0, oo, (40) 

m J J 0 


on equating the coefficients of k 1 , and as R 2 - R 2 = 2R 2 (h - h ) + o {fR 2 } . 
The quantities G. are obtained in this manipulation as 


G 0 = Ag 


(41) 
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X Q dcr, 


(42) 


R„ 


II 


h - h 


with more complex expressions for higher values of i (ibid, p. 122), being of the 
form 


G. = G i (h,h p .X 1 ,X 2 X ( i-i))- W 

Equation 38 reduces to Equation 40 when h = h p = 0 and 0 = 0. The Molodenskii 
problem is equivalent to Stokes' problem in such a case, the solution of which is 
Equation 15, which, on taking zero degree effects into account (Mather 1971c, 
p. 85) can be written as 


N f = 7 tv d + W 0 - U 0 ] = 


«o - U 0 

7 


R m M{Ag} R 


7 


4 tt y 


II 


f(0)Agdcr, (44) 


where M {A g} is the global mean value of Ag. 

The substitution of Equation 25 in its appropriately modified form, in Equation 40 
gives 


X: 


1 

217 


G + 


3V d 

2R 


(45) 


on adoption of the representation 


^ A f f , . x i 

V, = E - E k 77 

i =0 i=0 J J 


dcr + O { f V d ) 


(46) 


The second equality in Equation 46 would be consistent with Equation 25 only if 
there were no topography. If this inconsistency were removed (Molodenskii et al 
1962, p. 123), the final expression for the height anomaly would be a series in 
embedded in the form set out at 44 with some topographic correction terms 
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whose effects are purely local in character and need only be considered in areas 
of rugged topography, being functions of r j" 3 , the series being obtained when k = 1. 
In this case, 

+ 4Vy JJ f (S^) G dor + T - (47) 


DO 

£ G , m 

i -=0 

and T are the series of topographic correction terms whose form is given by 
Molodenskii (ibid). Moritz (1966, p. 91) has given alternative forms for G 15 and 
shows that if gravity anomalies are linearly correlated with elevation, Gj reduces 
to the terrain correction. Thus the combination of Equations 42 and 45 gives 
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(49) 


the second equality being based on the assumption of linear height correlation of 
gravity anomalies (ibid, p. 88). 

Notes :- 

(1) This technique will be practically effective only if the contributions of the 
higher G i are significantly smaller than those obtained for i = 0 and 1. The 
evaluation of any particular G i presupposes a knowledge of all X- (j < i), 
which in turn, are defined through Equation 45. 

The solution is therefore iterative, and as the series in x- is theoretically 
infinite, it is desirable that 


X, = oUO' 1 } x,.! (50) 

for efficient practical evaluation. As G 0 is the gravity anomaly, the first 
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iteration is the free air geoid, which contains over 90% of the power in the 
solution. It should therefore require only three iterations to obtain a solu- 
tion to order e 3 in h d , if Equation 50 were satisfied. 

(2) There would be little difficulty in meeting this criterion if the ratio (h - h )/ 
r 0 = o {10 -1 } . As oceanic regions comprising 70% of the Earth's surface 
area and non- mountainous regions make little or no significant contribution 
to topographical effects, the magnitude of the effects would be small if the 
above criterion were satisfied. All topography with grades in excess of 5° 
pose problems in this respect when they occur within a few km of the point 
of computation, distant zone effects being rapidly submerged by the r 0 _3 term. 
Also see section 3.5. 

(3) Serious embarrassment is caused when slopes exceed 1/4 v. Divergent, 
series are obtained, making an iterative approach unstable. Discussions on 
the problems of convergence are available in the literature (Moritz 1970; 
Moritz 1972; Krarup 1972). For a further discussion, see section 3.5. 

(4) The quantity G can have no first degree harmonic, as the solution of Stokes 
problem forbids the existence of such harmonics. Consequently, the refer- 
ence ellipsoid used for computing normal gravity is situated at the center of 
mass of the mass distribution needed to produce values of gravity at the 
surface of the Earth which would give rise to a gravity anomaly distribution 
equivalent to that of G. The writer is not aware of a detailed investigation 

of this problem but it is unlikely that the net effect would exceed o {5 x’lO cm}. 

(5) The extension of this technique to orders of accuracy greater than that of 
the flattening, is possible in theory. Such a solution could be obtained on 
including all effects of relevant magnitude in Equations 27 and 28, and on 
allowing for the existence of the atmosphere, noting that Stokes' integral is 
strictly valid only if there is no mass exterior to the physical surface. In 
addition it is necessary to take into account the Earth's ellipticity, especially 
when utilizing the orthogonal properties of surface harmonics. 

3.3 SOLUTION FROM DISCRETE VALUES 

This technique was originally proposed by Bjerhammar who summarizes the 

problem as follows (Bjerhammar 1964, p. 14). 

"A finite number of gravity data (gravity anomalies) is given for a non spherical 

surface, and it is required to find such a solution that the boundary values for the 

gravity data (gravity anomalies) are satisfied in all given points." 
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The Bjerhammar problem is differently posed to that of Molodenskii and a 
different approach is used for the representation of the gravity anomalies at the 
surface of the Earth. Working on the basis that the representation of the surface 
gravity field can only be in terms of samples taken at discrete points at the 
Earth's surface, Bjerhammar proposes the interpretation of such data in terms 
of a set of model anomalies Ag* on the surface of a sphere of radius R b which 
is less than or equal to the polar radius of the best fitting ellipsoid. The appro- 
priate requirement in Earth space is that any point on the Earth's surface lies 
exterior to the sphere of radius R b (the Bjerhammar sphere), whose center is 
collocated with the geocenter. 

The technique of solution can be summarized as follows. The surface of the 
Bjerhammar sphere is partitioned into a grid, each element of which has a sur- 
face area R b 2 dff, and is represented by the model gravity anomaly Ag*, assumed 
constant over the area. The disturbing potential V dp at any exterior point P, 
whose geocentric distance, as illustrated in Figure 5, is R p , is given by 


V 


dp 



2n + 1 

n ~ 1 



P n0 ( COS 0) dcr 


(51) 


under conditions applicable to Stokes' problem. Ag* obviously cannot have a 
first degree harmonic and the possibility of satisfying this condition in conjunc- 
tion with the geocentric collocation of the Bjerhammar sphere is subject to the 
same arguments as outlined in note 4 to section 3.2. 

The observational data is in the form of gravity anomalies Ag as determined at 
discrete points at the surface of the Earth. Using the fact that Poisson's integral 
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R b (R p 2 ~ Rb 2 ) 

4 77 



dcr, 


(52) 


applies without approximation to any function H which is harmonic exterior to 
the Bjerhammar sphere, it is possible to define gravity anomalies A g at all ex- 
terior points, if a surface distribution of the data set Ag* were available on the 
sphere. Alternately, if Ag are the gravity anomalies measured at the surface 
of the Earth, the equation defining the Ag in terms of the Ag* is obtained from 
Equation 52 as 
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(53) 
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where r is the distance between P i on the earth's surface, with geocentric 
distance R i5 and surface element da j on the Bjerhammar sphere. Equation 53, 
called a discrete integral equation by Bjerhammar (1968, p. 6), can be treated 
as a set of observation equations which can be solved by standard techniques for 
the elements Ag*. The technique is subject to certain practical difficulties when 
tested on models with heavy point masses between the sphere and the Earth's 
surface (ibid, p. 67). In such cases, Bjerhammar advocates the use of the dis- 
turbing potential rather than the gravity anomaly in Equation 53, presumably by 
recourse to an iterative procedure. The validity of the technique hinges on 
whether the resulting disturbing potential at points P , at the Earth's surface due 
to the Bjerhammar system is identical with that due to the Earth. For a sum- 
mary of the proof of this condition, see (Bjerhammar 1969, pp. 452-6). 

The instability of the inversion procedure due to the nature of gravitational 
attraction and its susceptibility to large masses locally (e.g. , mountainous re- 
gions) led Bjerhammar to suggest that the more stable disturbing potential V* 
on the Bjerhammar sphere be used as an intermediary in the solution on the 
following lines (ibid, p. 498 et seq.). 

The disturbing potential V di at on the Earth's surface is given by Equation 52 
as 
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Ag is obtained from Equation 13 on considering terms greater than o {f Ag}, as 
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Si 



3 ( r < pi 2 - R b 2 ) 2 
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K, 2 


V* dcr + o{f Ag} . 


(56) 

Equation 58 is simplified by differencing V* from the value V* 0 , which is the 
value of V* at the point on the Bjerhammar sphere corresponding to P. It can 
be shown (ibid, p. 499) that Equation 56 can be transformed to 
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(57) 


+ o { f A gp } , 


which is a generalized version of the Molodenskii inverse of Stokes' integral 

(Molodenskii et al. 1962, p. 50). 

Notes : 

(1) The use of this system would, at first glance appear to be a prohibitive 
task. This is not the case as the terms being integrated are scaled by r 4 , 
and hence only limited regions need be considered around each primary 
point at which evaluations are made. Details of test in the West Alps using 
a 5° x 5° area with a 15° x 15 c buffer zone, with basic subdivisions of 5' x 5', 
are given by Bjerhammar (ibid, p. 508), an iterative procedure being used 
to recover the Ag*. 

(2) The intellectual elegance of the method is enhanced by its ability to combine 

all manifestations of the Earth's gravitational field into a single solution 
entity. It must be added that this same end can be achieved by using the 
methods proposed by Krarup (1969), though the problems associated with 

: practical implementation have yet to be tackled in the case of high precision 
determinations. 
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(3) The factors which have to be taken into account to extend the solution to 
orders smaller than that of the flattening are similar to those outlined in 
section 3.2(5). 

(4) The solution, like that from Krarup's method, is unique for a given distri- 
bution of data. This, of course, does not mean that the answer obtained is 
correct to the order of accuracy with which the problem is formulated. 
Data requirements for solutions of the boundary value problem are dealt 
with in section 4. 

(5) For completeness, the solution should incorporate terms of zero degree as 
in section 3.2. 


3.4 SOLUTIONS FROM GREEN'S THIRD IDENTITY 

Considerable work has been done in this field (e.g., Arnold 1959; Koch 1965; 
Moritz 1965; Mather 1971c). The basic integral used is Green's third identity 
which is obtained by the application of Green's theorem to two scalars r - 1 and 
W which is harmonic in the volume V e exterior to a surface S. On combining 
the gravitational and rotational effects (e.g., Heiskanen and Moritz 1967, p. 15), 
the final expression obtained for the gravitational potential (W p ) at a point P on 
the surface S, on the assumption that all matter is contained within S and rotates 
with constant angular velocity <*>, is 



3 

where V = 3 -^- i, i being unit vectors along the axes x of a Cartesian co- 

i 1 

ordinate system, and r the distance of the elements of surface area dS and 
volume dV A interior to S, from P. A similar expression is obtained for the 
potential U p due to a gravitating reference ellipsoid which has the same rota- 
tional characteristics of the Earth, on considering the identical surface S, which 
is that of the Earth, when 


U P = Su-UV. Rl) ds - 2 ^ JJJidV,, (59) 

N in both Equations 58 and 59 being the unit normal vector to S at dS. 
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Both equations hold exactly if U and W are harmonic exterior to S. This condi- 
tion requires that no matter exists on either system exterior to the Earth's 
surface. The practical consequences, which are of significance when resolution 
approaching the order of the flattening is sought as the end result of computa- 
tions, are the following. 

1. The reference ellipsoid must always lie within S. As S coincides with 
the ocean surface over 70% of the Earth, the reference ellipsoid must 
be smaller than the ellipsoid which best fits the geoid by an amount 
greater than the largest negative geoid undulation (o { 10 2 m}). 

2. There should be no atmosphere exterior to S if Equation 60 is to hold 
to accuracies in excess of the order of 10 -2 V d . 

3. Both the reference ellipsoid and the Earth are assumed to rotate with 
the same constant angular velocity oj . Irregularities in the Earth's 
rotation have to be allowed for as corrections to observations in in- 
stances where such magnitudes are of significance. For details, see 
sections 2.2 and 4. 

The practice to date has been to treat atmospheric effects as those which should 
be modeled and allowed for as corrections to observations prior to use in com- 
putations (e.g. , IAG 1970, p. 18). In a recent solution Mather (1973, p. 28 et seq.) 
formulated a solution of the boundary value problem to o { e 3 h d ) by separating 
the gravitational effects of the atmosphere from those of the solid Earth and 
oceans. 

In conventional solutions, the disturbing potential V d is obtained on differencing 
Equations 58 and 59, when 

V dp = % - U P = i //( v d 5 • S 7 - 7 V • NV d ) dS. (60) 


This equation is not valid to orders smaller than o {10~ 2 V d > as it assumes the 
geopotential W to be harmonic outside S. A function which does satisfy Laplace's 
equation exterior to S is the potential W' due to the solid Earth and oceans, 
which is related to W by the relation 


W' = W - V , 

a ’ 


(61) 


where V a is the potential of the atmosphere, which is of order 10- 6 W, and more 
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significantly, V a = o {10~ 2 V d ) . As such, it is desirable to construct a theory 
which allows for its existence in the course of the derivation. 

The final solution using this technique can only be obtained by iteration, the 
number of iterations required as in the surface layer method, being a function 
of the accuracy sought. Favorable conditions for the adoption of an iterative 
procedure are the following, 

a. A significant amount (>90%) of the power should be generated in the 
first iteration. 

b. The iterative procedure should have the ability to converge to the 
correct result. 

c. The number of iterations necessary for achieving the desired degree 
of resolution should be as small as possible. 

When surface gravity is the sole source of information, the only procedure 
available for obtaining an adequate first approximation to the height anomaly h d 
is the use of Stokes' approach. Fundamental to this technique is the assumption 
that the disturbing potential is harmonic exterior to and on the surface S, and 
therefore can be expressed in the form 


where 





n / 1 , 


A 


I*. 

m “ 0 


and 


(62) 


(63) 


A - P (sin <p ) [C cos mX + S sin m A.] , (64) 

the last equation being the standard expression for a surface harmonic. 

The adoption of this model enables the combination of the effects of the dis- 
turbing potential V d and its vertical gradient 3V d /3hon using Equation 13, 
thereby transforming this formulation to a representation of the observed quan- 
tity, the gravity anomaly Ag. Details of the problems involved in obtaining a 
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solution of the boundary value problem to order e 3 h d are dealt with by Mather 
(1973, p. 31 et seq.). To preserve flexibility in the formulation of results to 
any required order of accuracy, it is desirable to retain physical relevance in 
the derivation by constructing the integral at 60 such that the disturbing poten- 
tial V d is replaced by the quantity 


= v d - v a 


( 65 ) 


where V a is the potential of the atmosphere. 

A generalized solution which did not consider either the existence of the atmo- 
sphere or the fact that the potential of the geoid was not known, the latter being 
disregarded after due consideration as a quantity which correctly cannot be de- 
termined from gravimetric methods alone (Molodenskii et al 1962, p. 104), was 
formulated by Molodenskii in 1945 (ibid, p. 93). A specific solution was given 
by Arnold (1959), which could be written as 


h 


dp 



A g "7b' a tan ,8 a ) f (0) d a 


( 66 ) 


+ 



2 77 y 


ff 1 

dh 

i-o 3 

(h p - h) + R sin 0 -jj: 


do- , 


where t, a are the components of the deflection of the vertical, f(0) is Stokes' 
function, tan /3 a being the components of the gradient of the ground slope in the 
north and east directions, 


where 


dh 
d r 


cos A^ tan 3 a , 


a; = 


a; 



(67) 


( 68 ) 


and r 0 is given by Equation 35. A revision of the derivation to the order of the 
flattening showed that (Mather 1971c, p. 85) 
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l dp 


Wp - Up 
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M {Ag} 
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4 77 7 
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Ag f (i p) dcr 


(69) 


R 2 

m 

2 77 7 


J fe h P " h ) + R m sin $ “f") “ “ y ^°- tan &]' da + o(f h d } , 


if 



o {f } , 


(70) 


► 

and Laplace's equation were satisfied to o {f V 2 V d } at all points exterior to, 
and on S. Equations (66) and (69) would be equivalent if 


1 

2 



7 ia tan Pa f 0/0 dcr 



<f a tan /3 a dcr . 


(71) 


The effect of the terms common to the kernels of the integrals on either side of 
the equality in Equation 71 can be expected to arise from only 30% of the surface 
area of the globe. Significant contributions to h d will be restricted to only about 
5% of the surface area, being about one order of magnitude smaller than A g if 
=*o {10 -4 } and tan f3 = o{10 -1 }. The high probability of correlation be- 
tween the signs of g a and /3 a in regions where the latter has significant magni- 
tude, indicate that this effect is likely to be always positive. As h d = o {10 2 m), 
it is realistic to estimate the effect of the above term as o {5 x 10 cm), the 
effect being consequential if regions of mountainous topography occur near the 
point of computation. 

The first term in the second integral at 69 converges much more quickly with 
increase in r Q and can be treated as a purely local effect. The advantage of the 
solution at 69 over that at 66 is the fact that all terms due to the interaction be- 
tween the ground slope and the slope of the equipotential can be treated as purely 
local effects, giving solutions where the neglected effects do not have magnitudes 
much in excess of the order of the flattening. Another advantage of the solution 
at 69 is its unambiguous definition in Earth space as the Stokesian term defined 
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a contribution with reference to an ellipsoid whose center is at the geocenter of 
an Earth which had no atmosphere. 

The generalization to order e 3 in h d (i.e., ±5 cm) deals not only with the effect 
of the topography, the interactions between the slopes of the topography and the 
equipotential surfaces of the Earth's gravitational field, as well as the atmo- 
sphere, but is also in keeping with the physical characteristics of the scalar 
potential (Mather 1973). It also establishes the nature of the relationship be- 
tween the Stokesian term and the indirect effect without limitations imposed by 
the simplistic approximations permitted by the adoption of a lower order of 
accuracy and identifies the anomaly to be used in Stokes integral. The geometry 
of the solution is also specified in Earth space, as the center of the reference 
ellipsoid is located at the center of mass G' of the solid Earth and oceans, whose 
co-ordinates X e . with respect to a geocentric Cartesian co-ordinate system are 
given by (ibid, p. 26) 


X ei 


_Ma 

M 




(72) 


where X ai are the co-ordinates of the center of mass of the model adopted for 
the Earth's atmosphere, M a and M e being the mass of the atmosphere and the 
solid Earth and oceans respectively. The final formulae obtained in the solution 
referred to are summarized below. 


*d P 


N fp + N cp , 


where the Stokesian term N f is given by 
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W 0 -U 0 _ M(Ag c } j 
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7 
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I. 


f 0/0 A g r . dcr, 


(73) 


(74) 


y p being the value of normal gravity at P' in Figure 1, R being the radius of the 
Brillouin sphere whose center is collocated with the center of mass of the solid 
Earth and oceans G' and contains the solid Earth and oceans. The gravity 
anomaly A g c is defined by 


A g c - A g t + A g 2 . 


(75) 
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where 

3V g V a 

A Si = A s + Th" + 2 R"’ 

m 


V a being the potential of the atmosphere, and 


( 76 ) 


if 


2 V, 


A 




S2 


R 


3h 


(77) 



3 2 A g 
3 h 2 


o (e 3 Ag} , 


where 


dR - R - a(l - f sin 2 4> c ) - h + o(f dR) , (78) 


3 Ag 

3 h 


r * 


tan <& c 

N f 

If BA 6.1 
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i 
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R 

m 

- 2 + o 

R 2 

m 


(79) 


and 


c^ - f + m - 3f sin 2 4> c + o (f 2 ) . 


(79) 


The angle \p is computed in calculations to the order of e 3 h d from the geocentric 
latitude 4> c and longitude \ as 

> b - cos -1 [sin sin <£ + cos c/> cos <f> cos d\] , (80) 

' cep c l* p 

where 


d\ = K - \ p , 


(81) 
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the value without subscript referring to do-, and those with subscript p to the 
point of computation P. 

The indirect effect N cp is given by 


N. 


cp 



+ 


1 

2 77 7 P 


r r 



3Vd 

3x a 


tan/S a + 


x a tan J3 a 






(82) 
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if 



3 2 Ag 
3h 2 


o (e 3 Ag} , 


r being the distance between dcr and P, R being the geocentric distance, given 
by 

R - a ( 1 _ f sin 2 <£ c ) + h + o(f 2 R} , (83) 

and h the ellipsoidal elevation. The other expressions which need definition are 

tan f3 a = - y<f a tan f3 a * N f tan + o{f 2 Ag} , (84) 

x a 1 

the two dimensional Cartesian co-ordinate system x a having the same significance 
as in Figure 2, 

R dh 

— tan 0 a = — (1 + c x ) sin «/» — , (85) 

r A r 
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where dh/ dr is defined by Equation 67, 


cos \ T ^ @ ^ 


1 AR /I \ 

2 sin 8 sin ^P ~ cos [~ 2 p + 8 I 

6 = - tan -1 + o{f 2 tan 0} 

„ o . 1 AR / 1 , 0 \ 

2 cos S sin ~ 2 p + sin I -y </> + 8 I 
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~[r - R p cos & + 8)]- 1, 


2 ir 

CA ' (1 +■ c ^ 2 " ’ 


= VT 
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+ o {f 2 } , 
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and 


r - 2R sin 


(93) 


Notes: 


(1) The adoption of an iterative procedure to solve Equations 73 to 93 cannot 
be avoided. See note (1) to section 3.2 for background. Mather (ibid, p. 49 
et seq.) has suggested an iterative procedure which requires three itera- 
tions of the expression for N cp to achieve an accuracy of 5 parts in 10 4 
(i.e., o {e 3 h d », while the Stokesian contribution need be evaluated only 
once, but in two stages. 

(2) The first two terms in Equation 74 are of zero degree and are meaningful 
only when the surface gravity field is sufficiently well defined to give rise 
to an adequate definition of the global mean value of surface gravity anoma- 
lies. Present day solutions, which are heavily dependent bn satellite deter- 
mined low degree harmonics of the Earth's gravitational field, tend to ignore 
the effect of these terms. For a further discussion of zero degree effects 
which would be of consequence in solutions based on adequate distributions 
of surface gravity alone, see section 5. 


(3) It could be construed that the conditions attached to Equations 77 and 82 are 
a limitation on the development outlined above. The relevant terms which 
are omitted have a net differential effect of 


1 

4 77 y 


yR f 0/0 - 



(dR) 2 


B 2 A g 

9h 2 


dcr 


on h d which would be negligible if the quantity B 2 Ag/ B h 2 had random error 
characteristics over areas larger than 10 4 km 2 with magnitudes of order 
10 - 8 mgal m~ 2 , which is only one of magnitude smaller than that of B 2 y/ B h 2 . 

(4) The components of the deflections of the vertical p a are computed on the 

principles outlined in Equations 16 to 24. A series of expressions which 

include effects with magnitudes of the order of the flattening are given by 
Mather (1971c, p. 86 et seq.). Such expressions should be adequate for 
evaluations of h d to o {e 3 h d > , but fail if accuracies of this type are required 

from the deflections themselves. Extension would principally require the 

use of an ellipsoidal co-ordinate system and a more careful evaluation of 
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some of the higher derivatives of characteristics of the gravitational field 
which have been assessed as having insignificant effects on the height 
anomaly. 


3.5 CONCLUSION 

The formulation of a solution for the boundary value problem at the physical 
surface of the Earth, in contrast to Stokes' problem where there is no topography 
exterior to the geoid, calls for the evaluation of "topographical" terms which 
arise as a consequence of 

(a) departures of the Earth's surface from a level surface; and 

(b) elevation of the point of computation above or below the surrounding 
topography. 

The first effect has contributions with long wavelength which, on present assess- 
ment of geoid determinations, should not have effects in excess of o {50 cm} 
unless rugged topography occurs in the vicinity of the point of computation. The 
second effect is a purely local one as it is scaled by the factor r~ 3 , as seen in 
Equations 42, 57 and 69. 

The limitation in theory, of the surface layer method is the heavy reliance it 
places on the convergence of the series in G i? defined at 43, in a mathematical 
sense. It may appear to be paradoxical in practical terms, that the slope of the 
topography of distant areas, which at least to a first order of approximation, 
are in isostatic compensation, can affect computations of the height anomaly. 
These terms occur because the gravity anomaly used in computations and de- 
fined by Equations 10 and 11, reflect the mass distribution of the Earth as it 
exists. If, on the other hand, a suggestion similar to that made by de Graaff 
Hunter (1958) calling for the smoothing of the Earth such that slopes in excess 
of 5° did not exist, were adopted, there would be little to choose between the 
methods outlined in this section for the solution of the boundary value problem. 

In such a case, all the iterative methods would require only 3 iterations to 
achieve 5 cm accuracy in h, . 

The conclusion that a model had to be adopted for the topography was also 
reached by Moritz (1972, p. 49) after a detailed study of the convergence of 
Molodenskii series. This approach has recently become anathema to physical 
geodesists (e.g., Molodenskii et al 1962, p. 118) as it involves making assump- 
tions about the density of material comprising the upper layers of the Earth's 
crust. In contrast, the quantities Ag, h, tan /3 and B* Ag/B bd must be con- 
sidered those which can be observed. In this sense, the solution described in 
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section 3.4 exhibits favorable convergence characteristics, as the series involved 
are not "open-ended", but controlled in magnitude, being terms in a rapidly con- 
vergent power series in the parameter f (= o {10~ 3 }). In practical terms, how- 
ever, the higher differential coefficients d 1 Ag/ Bh 5 - are unlikely to be determined 
and the adequacy of Equation 69 will depend largely on the magnitude and wave- 
length of the series (i!) -1 h 1 (^Ag/Bh 1 ). Cumulative magnitudes of o {±0.5 mgal} 
with wavelengths of 100 km or less can be considered to be acceptable for solu- 
tions to order e 3 h d ,as discussed in section 4. 

Reverting to the question of smoothening the topography in order that grades do 
not exceed 10 _1 , the problems which have to be resolved if such a procedure 
were adopted as everyday practice in physical geodesy are the following. 

(a) The principles underlying the transfer of mass and their associated 
consequences should be clearly defined. 

(b) The question of assigning a density for each element of transferred 
will have to be dealt with, as the resulting corrections to observed 
gravity, will depend on the model adopted for these masses. 

If such procedures were deemed to be necessary, it would be mandatory to adopt 
a model for the Earth with surface slopes less than 5°. The transfer of matter 
to achieve this goal will change both observed gravity as well as the location of 
the center of mass of the physical system. The exact numerical value of the 
corrections made will depend on the principles adopted for the mass transfer. 
Physical geodesists advocating this type of approach will have to face up to the 
philosophical problem of which elements of topography to flatten out or fill. It 
is most important that a single model be adopted in order that the geodetic com- 
munity is not subject to a confusing variety of results which are not in agree- 
ment, not because of significant factors, but merely as a consequence of the 
adopted smoothening procedure. The limited surface gravity data available at 
the present time continues to keep the above problem in the area of academic 
interest alone. It is one in which continued discussion is to be encouraged. 

The solution of Molodenskii's problem by means of analytical continuation 
(Moritz 1969; Marych 1969) has not been dealt with as it has been shown to be 
equivalent to the solution using the surface layer approach (Moritz 1971). 

Another method which may prove to have some benefits is the use of numerical 
integration techniques, on which published material is hard to come by. 


43 



4. PRACTICAL CONSIDERATIONS 


4.1 INTRODUCTION 

Practical considerations fall into two distinct categories. The first concerns 
the optimum sampling of data in order that the necessary precision can be 
achieved in numerical computations. The second is the extraction of the most 
probable results from whatever (inadequate) data is availably. The second falls 
beyond the scope of this review and is covered elsewhere in this Symposium. 

One exception is the use of Molodenskii and Cook truncation functions to obtain 
the maximum information from satellite determined gravity anomalies and local 
gravity fields. 

The problem can be summarized as follows. Over 90% of the power in h d comes 
from Stokes' integral. Many regions exist where dense local gravity fields exist, 
but where beyond some limiting angular distance 0 O , the available gravity data 
from the analysis of the orbital perturbations of near Earth satellites on com- 
bination with whatever surface gravity exists, can be represented as a set of 
surface harmonics, of the type given in Equations 63 and 64. The Free Air 
Geoid at Equation 74 can be written as (Molodenskii et al 1962, p. 147) 
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where the gravity anomaly A g to be used in Stokes' integral is expressed by the 
set of surface harmonics 
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Q n is Molodenskii's truncation function given by 
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P n0 (cos 0) being the Legendre zonal harmonic. Values of Q n for various values 
of 0 O are given by Molodenskii and his co-workers (ibid, p. 150) to n = 8, 
de Witte (1967) to n = 25, and Hagiwara (1972) ton = 18. The computational 
efficiency of this method over the use of surface quadratures techniques for 
distant zone effects, in the present era where distant zone fields are heavily 
dependent for quality on satellite data, is a factor of 70 (Ojengbede 1973, p. 32). 
hi practice, only a limited number (at present, up to degree and order 20) of 
such harmonics are available and a rounding off error will occur in the compu- 
tations, due to the existence of a residual in the power spectrum of gravity 
anomalies, on adopting the surface harmonic representation. Molodenskii uses 
an elegant technique to show that the use of harmonics to n = 8 with surface 
gravity representations up to 0 O = 23° results in errors less than 2 m, while 
extension of surface gravity coverage to 0 O =35° reduces the truncation error 
to less than 50 cm (Molodenskii et al 1962, p. 164). 

Similar considerations apply to the computation of the Vening Meinesz contribu- 
tion to the deflections of the vertical using Cook's truncation function (Cook 
1950, p. 377), the equation equivalent to 94 in this case being (de Witte 1967, 
p. 455) 
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A a being defined by Equation 19, while Cook's truncation function q n is given by 
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The relationship between the functions Q n and q n has been established by Hagiwara 
(1972, p. 461) who gives a proof of the equivalence of the developments by Molo- 
denskii and Cook. On using the same values of n and 0 O described in the previous 
paragraph, Molodenskii shows that the truncation errors in g a are less than 
1.1 arc sec and 0.2 arc sec respectively in the two cases given. 
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It can be concluded with confidence that the use of truncation functions is capable 
of giving a resolution equivalent to the best astro-geodetic results, as is borne 
out by determinations in Australia (Mather, Barlow and Fryer 1971, p. 19) and 
the tests carried out by Ojengbede (1973, p, 41). 


4.2 SAMPLING THE GRAVITY FIELD AT THE SURFACE OF THE EARTH 

The overwhelming majority of surface gravity data available at the present time 
has been surveyed for geophysical purposes motivated by regional considera- 
tions. Such information has to be carefully screened before being put to geodetic 
use. There are problems that arise in the establishment of the value of observed 
gravity itself. Until recently, most gravity determinations of quality were made 
by differential means using gravimeters. It is now possible to carry out an 
absolute determination of g with a resolution of ±50 /xgal using a transportable 
apparatus (Morelli et al 1971, p. 17) while resolution at the ±3 /xgal level has 
been reported by the apparatus at Sevres, France (Sakuma 1971). 

It is all important in the first instance that all values of observed gravity are 
correctly referred to the unified gravity standardization network defined by the 
"International Gravity Standardization Network 1971" (IGSN71) or an equivalent 
global control network in order that datum discrepancies may be minimized if 
not eliminated. 

The solution of the boundary value problem requires an evaluation of the gravity 
anomaly. This requires a knowledge of 

(a) the geodetic latitude of the gravity station to 0.04 arcsec (±3 cm) 
for an accuracy of 1 p-gal; and 

(b) the geopotential difference AW with respect to the geoid to ±0.003 kgal m 
for a resolution of 1 ^gal in the gravity anomaly Ag, in addition to the 
requirement stated earlier for values of observed gravity. This also 
calls for the definition of a datum for the geopotential differences on a 
global basis and to some desirable degree of resolution. 

The status at the present time is as follows. While IGSN71 is available, it is 
most unlikely that any of the large gravity data banks are reliably connected to 
this network in toto at the present time. Most values of normal gravity are 
computed from regional geodetic co-ordinates of gravity stations which are 
unlikely to differ from geocentric values by more than 10 arcsec. Thus all 
values of normal gravity computed in a continental area covered by one of the 
regional datums (usually up to 5% of the Earth's total surface area) are subject 
to systematic errors not exceeding ±1/4 mgal. 
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The lack of a global datum for geopotential cannot cause errors much in excess 
of 1/2 mgal if the datum for elevations were based on at least one year's tide 
gauge readings, and as there is no evidence available at present which indicates 
that the magnitude of stationary sea surface topography is much in excess of 
±2 m. While these magnitudes appear to be small, their effect on the evaluation 
of Stokes' integral is significant, being systematic in character. 

Present day geoid computations from surface gravity data are therefore limited 
in effectiveness as a consequence of irregularly distributed data which could be 
subject to systematic errors due to the effect of inadequately defined datums on 
the data set used in the computations. The existence of such effects cannot be 
tolerated when the data is required for the determination of the geoid with the 
highest possible precision in studies of sea surface topography, whose magnitude 
is unlikely to exceed 2-3 m. The term sea surface topography refers to depar- 
tures of the ocean surface from an equipotential surface of the Earth's gravita- 
tional field and is partially due to salinity, meteorological and tidal effects. The 
magnitude of the residual departures on allowing for these factors, and termed 
stationary effects, can only be estimated from manifestations along coastlines 
which have been obtained by comparing the results of geodetic levelling with tide 
gauge readings. Departures which cannot as yet be explained, have been re- 
ported in Australia (Hamon and Greig 1973) and the United States (Sturges 1972) 
with slopes approaching or in excess of 0.1 arcsec. On balancing existing satel- 
lite altimeter technology against the oceanographic requirements, it would ap- 
pear that a ±10 cm resolution in the determination of the geoid is a desirable 
goal for this purpose (Williamstown Report 1969, 3-2). 

The criteria governing the factors which constitute a "desirable" representation 
of the gravity field for the solution of the geodetic boundary value problem, is 
dependent on the requirements for the solution of Stokes' integral which, as dis- 
cussed earlier, provides over 90% of the power in the representation. This 
would apply to any of the techniques of solution described in section 3. The fol- 
lowing is a summary of a recent look at this problem (Mather 1973, p. 53 et seq.) 

A suitable form of Stokes' integral for quadratures evaluation is 


N/ cm > = K £ n, 2 £ Mij fO/r-p (100) 

i j 

where Ag. . is the value of the gravity anomaly representing a n.°x m 0 square, 

K = 1.58 xio- 2 , (101) 
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and /n ij = cos 0 cij or sin 0 Aj depending on whether a latitude-longitude or 
azimuth-distance system of co-ordinates is used. Equation 100 would be ade- 
quate if the subdivision of the basic n° x n P square into N ( = ^-j ) m° x m° 

squares (m < n . ) , where the k-th such square will be represented by the gravity 
anomaly A g k at an angular distance 0 k from the point of computation P such 
that 


Ag k = Ag + Cgk ; F(0 k ) = FO/O + c^, (102) 


A g and F ( 0) being given by 

1 n n 

A e = ; F <« =■ IT ZZ (103) 

k=l k=l 

and the use of these smaller subdivisions in the quadratures evaluation in lieu 
of the n.° x n.° and the appropriate area mean, did not reduce the quadrature 
error to below the desired order of accuracy (o {e }). This would happen if 

N 

£ c gk Sk = o{e} - ( 104 > 

k= 1 

implying no correlation whatever between variations in f (4>) and Ag over the 
n.° x n.° area. While the function F(0), given by 


F (> p) - f (0) sin \p , 


(105) 


has predictable variations, Ag defies accurate prediction except over very short 
distances and under carefully controlled conditions. As gravity has to be sampled 
at discrete points, the quadratures approach makes a representation procedure 
mandatory. Consequently, some finite element of surface area has to be repre- 
sented by a single observation. It is useful to bear in mind that 


(a) the global gravity standardization network available at present has a 
station accuracy of ±0.2 mgal (Morelli et al 1971, p.p, 6); 

(b) errors in gravimeter ties seldom exceed ±0.2 mgal if performed with 
adequate instruments and any sort of minimal care; and 
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(c) geopotential errors of o {3 kgal m} give rise to a o {1 mgal } error in 
the gravity anomaly. 

A precision of ±1 mgal in the gravity anomaly is relatively easy to obtain in 
areas where the regional geodetic level network is reasonably dense. The 
gravity anomaly also undergoes changes with position at different points in the 
basic square it is expected to represent. This penchant is characterized by a 
quantity introduced by de Graaff Hunter (1935), called the error of representa- 
tion E (A g mn > for an n° x m° square, which in the case of a fully represented 
square, is given by 


( E < A e >„„> 2 


L 

: r 1 


( A gj ~ A g) 2 

_ _ 


(106) 


A reliable value for E { A g nm ) is obtained from N evenly spaced values of Ag. 
covering the n° x m° square, A g being the mean value of the gravity anomaly, 
given by 

= £ L Ag i ■ (107) 

i= i 

Several estimates of this statistical characteristic of the gravity anomaly field 
at the surface of the Earth are available in the literature (e.g., ibid; Hirvonen 
1956; Molodenskii et all 1962, p. 172; Mather 1967, p. 131). Samples which are 
available at the present time from different parts of the globe, reflect the flatter 
continental areas. E { A g} n in such areas is a function of square size and, in 
general terms, can be expressed by the relations 


r ±c x Vn K° < n < 5° 

E{Ag} n = J (108) 

L ±C 2 n n < %° 


for an n° x n° square, where n is in degrees and E { Ag > n in mgal, when C t = 12 
and C 2 4 3 x 10. It can also be shown that E ( A g) n is a function of unsigned 
groundslope | {3 | , with magnitudes which can be as much as 5 times as great in 
very rugged mountainous areas especially when n is small. As such variations 
are not a function of elevation but of ground slope, it is estimated that about 
2-5% of the Earth's surface will require values of C x and C 2 which are 
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significantly greater than these given above, for an adequate representation of 
variations in the gravity anomaly. 

The n umb er of terms involved in the quadratures evaluation is a function of the 
accuracy desired in the computation. If the requirements of sea surface topog- 
raphy determinations (1 part in 10 4 ) were to be met, it would be necessary, for 
estimation purposes, to restrict square; sizes to those over which the contribu- 
tions of the terms containing the second differential coefficient of F(0) were 
held to o {e 3 h d >. The required number of summations is o {10 6 }. The study 
of the propagation of systematic and random error characteristics through Equa- 
tion 100, under these conditions, shows that an adequate representation of the 
surface gravity field which would enable the achievement of an accuracy of 
±10 cm in the final result would be one which had an E { Ag} value of ±3 mgal, 
if the data were not subject to systematic error in excess of ±50 figal. Such a 
representation is afforded by a 10 km grid in nonmountainous areas. While the 
estimation characteristics of gravitationally disturbed regions are covered by 
the above figures, which assume that oceanic fields will have a similar tendency 
to vary as continental data, regions characterized with larger ground slopes, 
have significantly greater values of E {Ag}. It would be necessary to reduce 
the size of the grid in such cases to retain E { Ag } at ±3 mgal. The use of 
smoothening techniques described in section 3.5 would of course reduce these 
values. It follows that present day techniques for establishing surface gravity 
anomalies are adequate for the determination of sea surface topography. It is 
interesting to note that the station spacing required on the above basis, is already 
available over large continental areas, like the United States, Canada and 
Australia, at the present time. 

The consequences of systematic errors in A g which hold the same sign over 
considerable extents, have significant effects on computed values of h d . A sys- 
tematic error e^ g which holds its magnitude over a n° x n° area but has random 
error characteristics over larger extents, is shown to have an effect e Ns on the 
computed value of h d given by (ibid, p. 65) 


e Ns " ±° < K " n e Ag> - 


(109) 


where K" = 10, for e Ns in cm, n in degrees and e^ g in mgal. If e^ s were held 
at ±5 cm, the estimate of the magnitude of the permissible systematic error e^, 
which is inversely proportional to its wavelength, varies from o {±5 mgal} when 
n = 0.1° to o {±0.1 mgal} when n = 5°. 

Likely sources of systematic error have been listed at the commencement of 
this subsection. The following conclusions can be drawn. 
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(1) IGSN71 would be an adequate gravity standardization net for sea sur- 
face topography studies only if the station density were 1 per 5 x 10 4 km 2 
and the errors of adjacent stations were not correlated at the 200 Atgal 
level. Neither of these conditions is likely to be satisfied. An adequate 
net would be afforded by stations at which absolute determinations had 
been carried out to ±50 /j. gal resolution, and a representation of 1 sta- 
tion per 10 5 6 km 2 . In the interim, it would be advisable that all gravity 
data should be subject to randomzation procedures at the level of the 
precision of the gravity standardization network, prior to use in solu- 
tions of the boundary value problem. 

(2) Gravity anomaly information on each geodetic datum should be cor- 
rected for changes in normal gravity due to the datum not being geo- 
centric (ibid, p. 16). 

(3) The term "geoid" which is assumed to be synonymous with both the 
global datum for elevations as well as the "undisturbed" free level of 
the sea, should be defined on the basis of models which afford resolu- 
tion with an accuracy of ±10 cm. 

A possible problem of some significance in the determination of sea surface 
topography and other high precision determinations of h d , is the existence of 
the sea surface topography itself with not insignificant amplitudes (e.g., 3-4 m) 
and substantial wavelengths. The evidence for the existence of such phenomena 
is widespread but based on purely coastal phenomena, as obtained from levelling- 
tide gauge comparisons. Extended studies of the sea surface using short pulse 
high resolution satellite altimeters should go a long way toward clarifying 
whether stationary sea surface topography is merely a coastal phenomenon, and 
if not, the dominant wavelengths with which it is prone to occur. The existence 
of stationary sea surface topography with 4000 km wavelengths and 2 m ampli- 
tudes would cause errors of o {±1 m} in h d . While this estimate is based on 
the largest estimate of the phenomenon presently available, the existence of 
such an effect will require an iteration of the determinations of h d . The rapidity 
with which these iterations converge are more a function of the wavelength of 
the stationary sea surface topography than of its amplitude. 


5. GRAVITY AND EARTH SPACE 
5.1 GRAVITY AND SCALE 

A problem which requires careful scrutiny is the possibility or otherwise of 
defining a scale for Earth space from gravity determinations at the surface of 
the Earth. It must be clearly emphasized that the ensuing development excludes 


51 



consideration of satellite data which constitutes the basis of low degree 
representations of the Earth's gravity field at the present time. The problem 
could be stated as follows. Given an adequate distribution of determinations of 
surface gravity, how are effects of zero degree h d0 in the global distribution of 
height anomalies to be interpreted. This effect can be written as 


W 0 -U 0 M{A gc } 

h d0 = Z R Z + N c0 + o{fh d0 }. (110) 


on considering Equations 73 and 74, N c0 being the contribution of zero degree 
by the indirect effect N c . W 0 is not known and it is common practice to assume 
the first term to be zero. The second and third terms will have finite magni- 
tudes which could be made equal to zero by changing the value of GM, and hpnce 
M { A g c }. This would, of course change the value of U 0 , implying an equivalent 
change in the estimate of W 0 if h d0 is to remain zero. h d0 could only be forced 
to take zero value if this is justified by some external condition. 

A more realistic procedure is to establish the numerical value of h d0 by 
analyzing the differences 


v i 


-Mi 


h n " h 0 i 


( 111 ) 


where h 0i is the ellipsoidal elevation based on the same reference ellipsoid as 
used in the gravimetric determination, but from independent observations, e.g. , 
geocentric satellite solutions. The adoption of a value for G M used in defining 
the system of reference will in turn define a value for W 0 , on extracting a zero 
degree residual from 111 and using it in Equation 110. 

Any "improvement" in the value of GM obtained from gravimetric determina-. 
tions is based on the assumption that the potential of the geoid W 0 is equal to 
that on the ellipsoid of reference U 0 . The geoid is a physical reality, being a 
manifestation of the mass distribution which gives the observed phenomena at 
the surface of the Earth, while U 0 is defined by the values chosen for each of 
the parameters a, GM, w and f defining the system of reference. It has been 
deduced that the term (W 0 - U 0 )/y is approximately 3 m if the ellipsoid were 
one of best fit to the geoid and GM were the best estimate' available for the 
Earth (Mather 1971c, p. 98), provided the free air anomaly had no zero degree 
harmonic. 

Thus any deductions which can be drawn about scale from gravimetric deter- 
minations alone are subject to ambiguity, if restricted to a single epoch, it 
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would be prudent to refrain from distorting the physical characteristic of the 
Earth represented by the value of G M, as determined with reliability, by inde- 
pendent means. The effect of zero degree deduced from comparisons described 
by Equation 111, can preferably be used in Equation 110 to deduce the true value 
of W Q . The effect of zero degree is therefore fully accounted for in Equation 74, 
and values of h d will no longer give rise to effects of zero degree in Equation 111 

A second effect of importance is the term of zero degree obtained on studying 
changes in observed gravity determined by the use of absolute techniques to 
resolutions approaching ±1 yu.gal, as determined on specially designed observing 
platforms, well distributed around the Earth, between successive epochs in time. 
Such changes can be interpreted as either reflecting an expansion of the Earth, 
as measured within the framework of the velocity of light, or else a change in 
the value of GM. For a discussion see (Mather 1972, p. 15). 


5.2 GRAVITY AND GEODETIC REFERENCE SYSTEMS 

The preceding development has assumed that the Earth has a fixed mass distri- 
bution subject to some periodic changes due to effects like tides. Such a descrip- 
tion would only be adequate if observations referred to a limited period of time, 
such as a few decades. There is considerable evidence which appears to point 
to the large scale redistribution of at least the masses constituting the Earth's 
crust over very long periods of time, with the attendant possibility of mass 
variations at greater depth depending on the nature of the mechanism which 
could give rise to such crustal motions. 

A possible consequence of such mass redistributions could be a motion of the 
Earth's center of mass (geocenter) with respect to the Earth's crust. The 
analysis of high precision determinations of absolute gravity at a well distributed 
net of observing platforms as described in the previous subsection, could provide 
a means of recovering the motion of the geocenter between epochs on analyzing 
the first degree harmonic of changes in absolute g (ibid 1972, p. 15). It should 
be pointed out that a problem in filtering out short period effects due to meteoro- 
logical causes has to be overcome before results of reliability are likely to be 
obtained. Fortunately, an estimate of the same effect can be obtained on studying 
changes in geocentric position of a global network of laser tracking stations, 
using dynamic techniques, to provide a verification of the effectiveness of the 
determination. 
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5.3 THE ROLE OF GRAVIMETRIC METHODS 

IN EARTH AND OCEAN PHYSICS 

Until recently, it was generally held that gravimetric methods if used with 
adequate data, provided the only non- controversial technique for computing 
ellipsoidal elevations with the same resolution as that available from geodetic 
levelling, thus completing the definition of geocentric position of points on the 
Earth’s surface in three dimensions. Position determination at the present 
time has not provided resolutions which can confidently be claimed to be better 
than 1 part in 10 6 . 

It is now clear that the most precise determination of geocentric position is 
required primarily for studies in Earth and ocean physics, rather than for any 
direct engineering or technological purpose. It would not be exaggeration to 
state that resolution to 1 part in 10 8 would be the aim of geodetic techniques 
being currently developed for such schemes. While there is no clear indication 
that surface methods, subject to restrictions imposed by atmospheric uncertain- 
ties, can be improved to meet these goals, extra-terrestrial techniques, like 
laser ranging to near Earth satellites and VLBI, promise that such goals may 
well be attained in the foreseeable future. There also is no reason to doubt at 
this stage, that transportable versions of these systems could not achieve this 
same degree of resolution. 

It would therefore appear that, with the passage of time, there would be less 
use of geodetic levelling and the systems of reference implicit in its concept, 
for use in Earth physics. The exception is of course the study of the instan- 
taneous geocentric position of the ocean surface and the interpretation of these 
results for the study of ocean circulation. The determination of the geoid with 
the highest possible precision are a necessary prerequisite for such studies. 
Gravity information will still have to be assembled and anomalies computed on 
the basis of elevations referred to an equipotential surface, the most convenient 
being the geoid. 

Three matters of significance which should be closely studied before undertaking 
the task of assembling an adequate gravity anomaly field for computation of 
geoid heights to ±10 cm, are the following. 

(a) The definition of a physical model to serve as a datum for elevations 
with an accuracy which is not more than a factor of 3 less the highest 
precision sought in the geoid solution. 

(b) Techniques to be used for minimizing the effect of gravity base station 
errors on geoid computations. 
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(c) The question of whether it is necessary to adopt a model for the 

"surface of measurement" and, if so, the nature of an acceptable model 
and the procedure to be adopted in converting measurements on the 
Earth's surface to equivalent quantities on the model. 
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